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ABSTRACT 

Non-local, time-dependent convection models have been used to explain the location 
of double-mode pulsations in Cepheids in the HR diagram as well as the existence and 
location of the red edge of the instability strip. These properties are highly sensitive 
to model parameters. We use 2D radiation hydro dynamical simulations with realistic 
microphysics and grey radiative-transfer to model a short period Cepheid. 

The simulations show that the strength of the convection zone varies significantly over 
the pulsation period and exhibits a phase shift relative to the variations in radius. We 
evaluate the convective flux and the work integral as predicted by the most common 
convection models. It turns out that over one pulsation cycle the model parameter <a c , 
has to be varied by up to a factor of beyond 2 to match the convective flux obtained 
from the simulations. To bring convective fluxes integrated over the He II convection 
zone and the overshoot zone below into agreement, this parameter has to be varied 
by a factor of up to ^ 7.5 (KuhfuB). 

We then present results on the energetics of the convection and overshoot zone by 
radially symmetric and fluctuating quantities. To successfully model this scenario by 
a static, one dimensional or even by a simple time-dependent model appears extremely 
challenging. We conclude that significant improvements are needed to make predictions 
based on ID models more robust and to improve the reliability of conclusions on the 
convection-pulsation coupling drawn from them. Multidimensional simulations can 
provide guidelines for developing descriptions of convection then applied in traditional 
ID modelling. 
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1 INTRODUCTION 


In the first paper of this series, |Mundprecht, Muthsam,~fe] 
Kupka ( 2013| (paper I henceforth), we have described exten¬ 
sions of the ANTARES code ( |Muthsam et al.||2010 ) which 
enable the investigation of the pulsation-convection inter¬ 
action by numerically solving the equations of radiation- 
hydrodynamics with realistic microphysics and opacities. We 
now turn to the description and analysis of a 2D nonlinear 
Cepheid model. While the first paper centered around tech¬ 
nical issues, the present one concentrates on the physical 
properties of the convection zone caused by the ionisation 
of He II which also hosts much of the pulsational driving ex¬ 
hibited by the ^-mechanism. To motivate this work we first 
discuss here a few items related to the physics of Cepheid 
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pulsations which explains the interest in the modelling of the 
pulsation-convection interaction. Rewards from such work 
are expected in the course of time, both in the immediate 
and in the later future. 


The success of the effort put into modelling the 
pulsation-convection interaction has, to a large degree, been 
measured by the accuracy up to which properties of indi¬ 
vidual Cepheids can be reproduced (light curve, in partic¬ 
ular also special features such as humps and bumps, etc.). 
There is now a large body of ID simulations which are uti¬ 
lized in various astrophysical contexts. Some problems and 
trends in Cepheid research have been described in two re¬ 
views, iBuchlerl (119971) and iBuchlerl (|2009|. The status and 


progress in research on the related RR Lyr stars has recently 
been reviewed by Marconi (2009). It is remarkable that both 
authors stress open problems in the treatment of pulsation- 
convection interaction. Instead of repeating the problems 
discussed in those reviews we want to shortly describe a 
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few problem areas which are either different from those dis¬ 
cussed there or are considered here from a different point of 
view, and where better insights to the pulsation-convection 
interaction can expected to be crucial for progress. 

Let us, first of all, point out that in the usual ID sim¬ 
ulations for those stars convection is modelled by simplified 
non-local, time-dependent models, specifically adapted to 
an environment varying due to the radial pulsations. Gen¬ 
erally, such models need some equation for the turbulent 
kinetic energy, and due to the variable environment at least 
a time evolution equation for the turbulent kinetic energy is 
necessary. This is in contrast to the possibly static equation 
for the turbulent energy which is applicable in stars which 
do not exhibit pulsation. See the discussion in |Buchler X\ 


Kollath (2000) where the most widely used models are de¬ 


scribed. 

These models require, however, to know a number of 
parameters, the so-called cds, which are unknown a priori. 
Typically, these parameters are calibrated so that models re¬ 
produce properties which those objects exhibit, for example, 
the existence and the location of the red edge of the pulsa¬ 
tion instability strip. Sometimes they are simply guessed. 
Since the number of these parameters is, in addition, con¬ 
siderable, and even the form of the added equations or clo¬ 
sure approximations may well be ambiguous concern must 
be expressed as to the reliability of this approach. 

As a consequence, these models must be tested and, if 
feasible, the parameters calibrated in different, independent 
ways. Presently, it is not known to what extent the results 
obtained from conventional models are affected, at least in 
the case of cooler pulsating stars where effects due to convec¬ 
tion are largest. It is therefore of interest to consider a spe¬ 
cific problem where what had once been considered progress 
has again been challenged just on grounds of inappropriate 
modelling of convection: the (ID) simulation of double mode 
Cepheids. 

In parameter regions where double mode behaviour is 
observed, nonlinear, radially symmetric models for pulsating 
stars have failed to reproduce such a behaviour for a long 
time. It has therefore been considered to be an important 
breakthrough when, ultimately, models did exhibit double 
mode pulsation indeed. These results have, however, been 
criticized by Smolec X Moskalik (2008), where also refer¬ 
ences to other papers are given, on grounds of inappropri¬ 
ate modelling of convection, in particular because of much 
too large overshoot. Indeed, the question of how to explain 
double mode pulsation in Cepheids is now considered open 
again. 

Quite generally, these closure parameters of the convec¬ 
tion models are known to depend on effective temperature, 
chemical composition, etc., insofar as they are can be deter¬ 
mined at all with some reasonable reliability. This holds al¬ 
ready true for the ” classical” mixing-length parameter l/H p . 
Consider for example the discussion in |Bono, Castellani, 


Marconi (2002). Increasingly, Cepheids and related variables 


are used to investigate stellar populations, see, for example, 


Marconi 


(2012). Hence, a proper assessment of the validity of 


conclusions is becoming both more important and, because 
of the possible influence of differing chemical abundances on 
the mixing length and other parameters, more difficult. 

Another emerging area of research is the detection of 
multiple periodicities present in some of these objects due 


to the advent of long, uninterrupted observations of high 
quality with the CoRoT and Kepler missions. The detection 
of nonradial modes, e.g. Guggenberger et al. (2012), amounts 
to new tasks for nonlinear modelling and in particular leads 
to the need for multidimensional simulations. 

Let us also point out the usefulness of multidimensional 
Cepheid and RR Lyr studies for investigation of the atmo¬ 
spheres of these objects and the implications of such work. 
Whereas in solar and in general stellar physics use is in¬ 
creasingly made of multidimensional model atmospheres for 
abundance analysis and other types of research, e.g. |Asplund| 


et al. (2009), static model atmospheres are predominantly 


used in work about Cepheids with only few exceptions, e.g. 


Nardetto et al. (2004). An appropriate treatment of the so- 


called projection factor which relates actual pulsation ve¬ 
locities to those observed turns out to be of importance in 
establishing the rotation law of the galaxy and issues regard¬ 
ing Cepheid distance scale calibration (cf. |Nardetto et al. 


2009). Quite generally, the increasingly sophisticated analy¬ 


sis of light curves of Cepheids and RR Lyr stars, e.g. |Kanbur| 


et al. (2010), should be backed up by a comparable sophis¬ 


tication of atmospheric models, therefore possibly hydrody- 
namical and, if convection is important, multidimensional 
models. 

With these long-term goals in mind we present in this 
paper 2D hydrodynamical models of a Cepheid. As discussed 
already in paper I, the vastly different spatial scales and 
therefore stringent resolution requirements make this a com¬ 
putationally quite intensive task even in 2D. In the following 
we investigate properties of the He II partial ionisation zone, 
in particular issues regarding the energetics of the convec¬ 
tion and the overshoot zone and implications for modelling 
of convection done in the traditional way in this area (mod¬ 
els of such a sort that they fit into the framework of the 
geometrically one dimensional hydrodynamic modelling of 
radially pulsating stars), and we discuss issues of the atmo¬ 
spheric structure and numerical complications which are not 
so severe in low resolution and show their full relevance only 
in high resolution simulations. 

The remainder of this paper is organised as follows. In 
Sect. [2] we describe the physical parameters of the Cepheid 
which we model with our simulations. The analysis of the en¬ 
thalpy flux in and around the partial ionisation zone of He II 
and the mechanical work performed by the pulsation over a 
complete cycle as well as the comparison of these quantities 
with predictions from time-dependent convection models is 
presented in Sect. [3] In Sect. [4] we provide a discussion of 
our results and present our conclusions. 


2 MODELS 

The model we discuss in this paper corresponds to a fun¬ 
damental tone Cepheid. Its parameters have already been 
given in paper I. We repeat them here for convenience. 

Our star has an effective temperature of T e s = 5125 K 
and a gravitational acceleration at its surface of log g ~ 1.97. 
Its luminosity is L ^ 913 L®, the initial radius is R — 
26.8 Gm (which is about 38.5 R©), and the total stellar 
mass amounts to 5 M®. 

The chemical composition is characterized by X — 0.7 
and Z — 0.01, the mix for Z being taken from Crevesse X\ 
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Noels (1993). The outer 11.3 Gm are modelled. As a result 


the pulsation period P measured from our simulations is 
about 3.85 days which is slightly shorter than if one included 
the entire star in the computational domain. 

The present work is based on a model whose sector has 
an opening angle of 10°. The vertical cell size ranges from 
0.47 Mm at the top to 124 Mm at the bottom. Further de¬ 
tails on the numerical parameters for the model and on the 
numerical methods used to solve the dynamical equations as 
well as a description of the dynamical equations themselves 
are given in paper I. 


3 PULSATION AND CONVECTION 

In this section we discuss properties of the He II of the model 
just mentioned. This work is based on data covering ten 
useful pulsation periods after model relaxation. Averaging 
over ten periods allows for a sufficient removal of special 
properties of individual periods. After all, a problem typi¬ 
cally experienced when running 2D convection simulations 
the small number of plumes present at any given time un¬ 
less the domain is very broad. This makes for poor statistics 
compared to 3D runs with similar horizontal extent. The 
He II zone is well resolved in each of the models we consider 
in the following. 


3.1 Comparison to traditional convection models 

3.1.1 Definition of the fluxes 


Given the importance of convection modelling in traditional 
(ID) Cepheid calculations naturally great efforts have been 
undertaken to derive such models. The most frequently used 
models are those of KuhfuB (1986) and of Stellingwerf (1982) 
as well as a number of variants based on them. 

Convection modelling in the sense of KuhfuB has 
been adopted with modifications in the code described by 


Wuchterl & Feuchtinger (1998). Variants of the Stelling¬ 


werf model and the KuhfuB model have been used in 


codes such as those of Bono & Stellingwerf| (1994) and 


the Florida-Budapest cooperation (e.g Buchler, Kollath, & 
Marom (1997)) and also in the more recent code described 


Smolec &; Moskalik (2008). 


These convection models introduce one additional time- 
dependent differential equation for the turbulent kinetic en¬ 
ergy density, e t , plus algebraic recipes for various physi¬ 
cal quantities. More elaborate convection models using the 
Reynolds stress approach have been developed, in particular 
those of Canuto (1997) and Canuto & Dubovikov (1998), but 
they have not yet been integrated into ID nonlinear stellar 
pulsation codes. The comparably complex model of |Xiong| 


& Deng (1997) has been used to study the pulsational sta¬ 


bility of red giants, but only at the level of providing the 
background model for linear, nonadiabatic stability analy¬ 
ses (see jXiong &; Deng (1998) and Xiong &, Deng (2007)). 
We therefore concentrate on examining the prescriptions of 
the KuhfuB and Stellingwerf models. 

One of the differences between these models is the way 
in which the convective flux is expressed from the basic 
quantities, i.e., those derived from the Navier-Stokes equa¬ 
tions plus e t . Errors in the prediction of the convective flux 


may stem from two sources: from errors in et itself (and 
therefore the basic equation for this quantity) and from an 
inadequate algebraic expression used to calculate the con¬ 
vective flux F c . In the tests performed here we compute et 
from our 2D simulations and thus we test the algebraic form 
for F c . For two reasons this is the proper way to proceed in 
the present context. The first is a physical one and the other 
is of a primarily numerical nature. 

The physical reason is conected with the formulation of 
the time dependent convection (TDC) equation. In the form 
given by Buchler &; Kollath ( |2000| it reads 


dtet + (pt+pAdtV = -Prd r (r 2 F t ) + C, (1) 

pr z 

where pt is the turbulent pressure and p u is the viscous eddy 
pressure. F t is the flux of turbulent kinetic energy and C is 
a coupling term of the form 


C = S-e , 


( 2 ) 


where S denotes a source term and e is the dissipation rate of 
the (turbulent) kinetic energy. Both the viscous eddy pres¬ 
sure p u and the dissipation rate e concern actions on a very 
small scale and in actual practice cannot be computed with¬ 
out model assumptions concerning just those small scale 
processes. From the standpoint of 2D or 3D modelling that 
means, however, that these quantities cannot readily or un¬ 
ambiguously be recovered from the numerical simulations. 
After all, the latter contain information about the resolved 
scales only. The hope and expectation (which to a consider¬ 
able extent was supported e.g. by simulations of solar gran¬ 
ulation) in connection with such models is that they deliver 
useful and reliable answers even if the parameters (Reynolds 
number, Rayleigh number, etc.) in the numerical models dif¬ 
fer widely from those found in actual objects (cf. [Kupka 
(2009)). At the same time it should be kept in mind that 
Eq. § in itself is phenomenological only. It is in at least 
one way connected with simplified expressions for e and con¬ 
nected to the Kolmogorov dissipation picture as described 
Buchler & Kollath (2000). However, Eq. |l]) is not built 


into the multidimensional simulations. 

Considered from a more basic point of view, Eq. 0 is 
not conservative. Moreover it cannot be brought into that 
form with the approximate associated expressions for e and 
S . The physical reason for this is that turbulent kinetic en¬ 
ergy considered on its own is simply not a conserved quan¬ 
tity: there has to be a source of turbulent kinetic energy 
to prevent it from decaying. In the models used in practice 
the source and dissipation terms are, from a quite funda¬ 
mental point of view, not provided in conservation form. On 
the other hand, our numerical methodology in ANTARES 
rests on conservation form. This is for the very good reason 
that in hydrodynamics only conservation form can be ex¬ 
pected to deliver correct answers for shock speeds and the 
like (cf. [Trangenstein (2009)). Nevertheless, even if we chose 
to solve Eq. (|lj alongside, it would not have the same, or 
any, backreaction on our simulations because it is not one of 
the dynamic equations we solve. It would therefore concern 
a quantity different from that one appearing in the TDC 
models despite of a formal, superficial analogy. 


Therefore, similar to Gastine & Dintrans (2011), we un¬ 


dertake a direct comparison of the convective flux between 
the numerical hydrodynamical simulations and its approx- 
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imation provided by the TDC models. Each of the TDC 
based fluxes has the form 


F cs * = a c ,*r]*(j)* 


( 3 ) 


Here, * is KF for KuhfuB and SW for Stellingwerf. 77 * is a 
correction factor whose purpose is to account for radiative 
energy losses of the convective elements. For the KuhfuB 
and the Stellingwerf case a slightly different functional form 
is used (see below). 77 * contains, among others, the Peclet 
number Pe*, which, again, is calculated somewhat differently 
in the KuhfuB and in the Stellingwerf setting. 

Our use of a c complies with Gastine & Dintrans (2011). 


Note that e.g. Buchler &; Kollath (2000) use a somewhat 


different way to write the factor. 

The functions 0*, <a c ,* and 77 * depend on the depth. In 
the Stellingwerf model 0sw is computed from 


0SW — O • p • e t • sign(Vsupad) ’ l^supadl} 
where © = yj c p • T/V a d- It furthermore adopts 

1 


V sw — 


1 + <a r /Pesw ’ 


( 4 ) 


( 5 ) 


with a Peclet number computed from the following expres¬ 
sion: 


Pesw — H p - y/et 


Cp • n ross " P 
16(7 rp 3 

3 ' 1 


( 6 ) 


This form was proposed in Buchler & Kollath (2000) (see 
Sect. 5.3 therein) originally with the model of KuhfuB in 
mind, although at this level the expressions for the convec¬ 
tive flux have the same general structure .<© for both the 
Stellingwerf and the KuhfuB models. In © , a r is a param¬ 
eter, unknown from the outset, /^ roS s is the Rosseland mean 
opacity, a the Stefan-Boltzmann constant and T the tem¬ 
perature. Moreover, c p is the usual heat capacity at constant 
pressure, p denotes the density, e t the turbulent kinetic en¬ 
ergy density, and H p stands for the pressure scale height. 
Finally, V sup ad is the superadiabatic temperature gradient, 


^supad — ^ ^ad? 


( 7 ) 


and V denotes the actual temperature gradient d In T/d In P 
while V a d is the adiabatic temperature gradient. 

For the Kuhfufi model we have 


0KF — Cp - p y/~Ct ' Vsupad- 


( 8 ) 


To demonstrate the effect of radiative losses on the KuhfuB 


model if approximated as in Buchler & Kollath (2000) we 


introduce here a factor 77 KF which contains some parameter 
a T . Again, it influences the efficiency with which the model 
diffuses surplus energy of convective elements and it is as¬ 
sumed to be of the form 

OLv • PeKF , n x 

m? = tz -si—- ( 9 ) 

1 + <a r • PeKF 

The Peclet number for the KuhfuB model is computed from 


\/0A -tv Cp- p 

PeKF = --- • Hp - yfef - ——, 

O A ra d 


( 10 ) 


where Arad = ( 16 crT 3 )/( 3 ft r 0 ssp) is the radiative conduc¬ 
tivity. This is the expression for Pe derived by Can uto k\ 
Dubovikov ( 1998 ) for the case that the dissipation rate e is 


l — Hp. Note that Eq. § can be transformed into Eq. 
the role of a r is reversed between the two, i.e. 


dr,KF — QJ r) SW 


5/(V04 7 


( 11 ) 


while PeKF = (VWAtt/5) Pesw- Consequently, all the dif- 
ferences between Eq. <© and Eq. § can be absorbed into 
the numerical value for a r and the two forms of comput¬ 
ing 77 are equivalent. Eq. © has the advantage of a more 
straightforward comparison with the model of |Canuto ~fe| 
Dubovikov (1998) and extensions thereof in future work. We 


note that more detailed prescriptions to compute radiative 
losses within the KuhfuB model itself have been suggested 


by Wuchterl & Feuchtinger (1998), which, however, quite 


closely couple the computation of the convective flux with 
the dynamical equation for et. In this study we limit our¬ 
selves to the prescription given by Eq. © which can readily 
be applied to both TDC models, as intended by [Buchler fe| 


Kollath 


( 2000 ). 


computed by means of a mixing length l in the special case 


Last but not least we also compute the convective flux 
from our simulations through 

F c =u' r -(p-h). (12) 

Here u' r = u r — up is the vertical convective velocity (the 
radial velocity after subtraction of the radial pulsational ve¬ 
locity) and h the enthalpy. 

3.1.2 Flux properties and comparison 

In order to compare the fluxes used in traditional Cepheid 
modelling with those from our simulations and relating them 
to the actual convective fluxes as calculated form the sim¬ 
ulations we have averaged the latter ones over ten useful 
pulsation periods. In this way we improve the statistics. Af¬ 
ter all, it is known that 2D simulations have disadvantages 
regarding statistics compared to 3D simulations due to the 
small number of convective plumes present at any given in¬ 
stance of time unless the domain is very wide. For each of 
the results we discuss in this section the fluxes have been 
averaged horizontally. 

Fig. ^displays the convective flux plotted against phase 
(averaged over ten periods). The variation of convection in 
intensity as a function of the phase is clearly visible. As 
Fig. 1(a) shows, convective activity (as witnessed by con¬ 
vective flux) commences its activity high up in the convec- 
tively unstable layer. Subsequently, plumes descend down¬ 
wards. While their effect is still recognizable deep down (in 
the overshoot zone, in fact), new convective activity again 
sets in at the top of the He 11 partial ionisation zone. 

The two subfigures in Fig. ^ display two variants of 
the convective flux, namely F c3 — u' r - (ph — ph ), based on 
entropy perturbation again st the mean, in Fig. |l(a)[ and 
F c — u r - (ph) in Fig. l(b)| where the total enthalpy enters 
(cf. §. In the upper picture there is a noticeable orange- 
coloured streak, corresponding to plumes directed down¬ 
wards, in particular in the overshoot region. They have a 
large negative value of velocity, and ph < ph because they 
are hotter than the ambience in the overshoot region. This 
leads to a large positive value of F c 3 — u r - (ph — ph ), but 
to a negative value for F c . 

In Fig.[2]the mean convective flux F c is compared to its 
approximation provided by the TDC models over the total 
computational domain. The left panel refers to the KuhfuB 
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0 0.5 1 


phase 

(b) F c =u' r - ( ph) 


F c3 

5,8e+06 

■4e+06 

2e+06 

-l,2e+06 


F c 

^,82e+06 

^2,5e+06 


^-2,5e+06 

-2,54e+06 


Figure 1. The convective flux obtained from the simulation (in colour) after averaging over ten periods is plotted here against phase. 
The vertical axis shows the radial extent of the model, while phase is plotted horizontally. In the upper picture the convective flux F c 3 
is computed only from the enthalpy perturbation. The units are Js - 1 m -2 for both pictures. Note the asymmetry of the radius as a 
function of phase, a feature well known through observations of Cepheids. Convective entities form high in the He II partial ionisation 
zone are found to migrate downwards into the overshooting region. 


and Stellingwerf fluxes without the correction factor, which 
amounts to 77 * = 1. In the panel to the right the correction 
factor is included. Similiar to the findings in Gastine & Din- 
trans (2011), but less pronounced, the KuhfuB formulation 


overestimates the convective flux in the overshooting region 
also in our case. Partial remedy can be achieved by including 
the Peclet-factor 77 * as is evident from Fig. |2(b)| The results 
based on the Stellingwerf model turn out to be only weakly 
dependent on the Peclet correction. - The approximations 
also predict a slight negative overshooting above the con¬ 
vection zone. The transition from convection to overshoot¬ 
ing occurs closer to the surface than for the convective flux 
proper. - In all figures and tables hereafter the Peclet-factor 
will therefore be included. 


Fig. [3] shows a c ,* as a function of phase for one full pe¬ 
riod (actually even from phase - 0.2 to 1.2 for clearer repre¬ 
sentation). a c ,* is obtained from a c ,* — F c ,*/ (4>* • rj*). The 
reddish vertical stripes represent the upper and lower limits, 
respectively, of the He 11 convection zone proper quite well. 
A systematic variation of <a c ,* with both the depth in the 
convection zone and with pulsation phase is immediately ap¬ 
parent. Furthermore, comparing corresponding values from 


the convection zone and the overshoot zone, one can easily 
spot phases where the <a c ,*-values differ by a factor 5-15 and 
hence by an order of magnitude. 

We therefore need to investigate the coefficients a c ,* 
as a function of phase separately averaged over the convec¬ 
tion zone and the overshooting zone. To this end the follow¬ 
ing procedure has proved reliable and useful to determine 
whether a depth point belongs, at a given instant of time, 
either to the convection zone or the overshooting zone or to 
neither of them. 

A layer belongs to the convection zone if, at this depth, 
more than 90 percent of the horizontal points have a pos¬ 
itive superadiabatic gradient, V — V a d > 0. The gradients 
are computed from local values. Points below the convec¬ 
tion zone belong to the overshooting region if they satisfy 
r 3 < 1.51 and fulfill, in addition, V — V a d < —0.03. This 
choice guarantees both for the convection zone and the over¬ 
shooting zone that the upper and lower boundary of the con¬ 
vection zone are excluded. Near these boundaries the con¬ 
vective flux changes sign and hence is practically 0 so that 
the ratio F c */ ((ft* • 77 *) requested for the evaluation of a c 
(Cf. Fig. [3f leads to meaningless results. Moreover, in this 
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250 300 350 400 450 



(a) Comparision of F c and </>* 


250 300 350 400 450 



Figure 2. Mean convective flux in the simulation compared to the flux values for the Kuhfufl and Stellingwerf expressions for the 
convective flux. 77 = 1 is assumed in the left panel whereas 77 is determined according to Eq. {H and Eq. § in the right panel. The top 
of the simulation box (at gridpoint 1 , not in the range of the figure) is to the left of each panel. 




phase 


(a) KF - convection zone, a r = 0.3 


0 0.5 1 



0 0.5 1 


phase 


(b) SW - convection zone, a r =0.1 


QSW P a^t,2 



Figure 3. Horizontally averaged a c -values for the Kuhfufl (left panel) and Stellingwerf (right panel) expressions for the convective flux 
including 77 * and including the y-component of the turbulent kinetic energy (e^), indicated by a colour scale from blue (small values) 
to red (large values). The horizontal axis displays phase while depth is plotted along the vertical axis (range of radius from 24.5 Gm to 
26.4 Gm). An average over ten periods is shown. 


region the enthalpy flux is small and the ratio F c ,*/ ( 0 * • 77 *) 
is physically unimportant. Inspection of the simulation data 
shows that the overshooting region defined in this way nicely 
contains the downwards penetrating plumes. 


3.1.3 Domains of convection 

Given the fundamentally different nature of the convection 
zone proper and the overshooting region it is mandatory to 
investigate from the outset the achievements of the convec¬ 
tion models for those two zones separately. To obtain a basis 
for those comparisons we have averaged data, which them¬ 


selves are horizontal averages at each instant of time, over 
ten periods as already explained in Sect. [3TT72| 

The comparison of the convective flux obtained from 
the simulations with the TDC prescriptions calculated as de¬ 
scribed in Sect. |3.1.11 is summarised in Fig. [4] and in Table [I] 
The meaning of the abbreviations designating the variants 
of the convection models are given in the caption of Table ^ 
and a c is determined via Eq. §• The large differences be¬ 
tween the convection zone and the overshooting zone clearly 
stand out as does the strong dependency of the values for 
a c on the adopted value for a r . In addition, there is a strong 
variation of a c with phase. The dependence of a c on a r is of 
opposite sense in the Kuhfufl and the Stellingwerf model, re- 
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o 0.5 1 



phase 

(a) KF - convection zone 


0 0.5 1 



(b) SW - convection zone 


o 0.5 1 



(c) KF - overshooting region 


0 0.5 1 



(d) SW - overshooting region 


Figure 4. The parameter a c as a function of time: horizontal averages over ~1 period. In the different panels KF denotes results from 
testing the KuhfuB model and SW results from testing the Stellingwerf model. Letters r, s, a, and b after P in the designation identify 
different values of a r . The results denoted by et,i include only the vertical velocity in the computation of the turbulent energy while for 
results denoted by et ,2 the horizontal velocity was included as well. 


TDC 

variant 


convection zone 

overshooting region 


designation 

Or 

a c 


a max 

«c 


a max 


O 

a min 

o 

a min 


Pr e t ,i 

0.08 

11.2 

31 

2.8 

1.59 

10 

1.4 

KF 

Pr e t ,2 

0.08 

7.2 

22 

2.3 

0.95 

14 

1.7 

Ps e t ,i 

0.30 

3.2 

30 

2.7 

0.46 

10 

1.4 


Ps e t ,2 

0.30 

2.1 

22 

2.2 

0.28 

14 

1.7 


Pa 6t,l 

0.10 

1.13 

26 

2.2 

0.48 

7 

1.3 

SW 

Pa et ,2 

0.10 

0.74 

15 

1.5 

0.26 

10 

1.4 

Pb e t,l 

0.25 

1.33 

28 

2.3 

0.55 

8 

1.3 


Pb e t,2 

0.25 

0.84 

16 

1.6 

0.30 

10 

1.4 


Table 1. a-values for various convection models: the designation indicates whether the ^-component only 
is used for evaluating turbulent kinetic energy (et,i) or whether also the y-component is used (et,2)- The 
letters r, s, a, and b in P r , P s , etc. refer to the value of a r adopted (listed in the next column). a c is the 
mean value, a the standard deviation in percent of the mean value (averaged each over ten periods). The 
ratio Q:max/Q: m in is a measure for the total amplitude of the variation of a c over one period. 


spectively. This is an immediate consequence of the different 


definitions of a v in both cases (see Sect. 3.1.1). 


The quality, i.e. the constancy with phase, of a certain 
set of cFs may be judged either from their standard devia¬ 
tion a or from the ratio a m ax/amin (Table]!]). In all cases the 
two criteria lead to the same assessment of quality. In addi¬ 


tion, there is the possibility to either calculate the turbulent 
kinetic energy from the vertical velocity alone or from both 
velocity components. While the second option may appear 
to be the more natural one, the first option is of interest be¬ 
cause some authors (for instance Gastine & Din trans|[2011 ~ 


specifically Eq. (2) in their paper) consider the contribu- 
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0 0.5 1 



phase 

(a) KF - convection zone 


0 0.5 1 



phase 

(b) SW - convection zone 



Figure 5. The parameter (3 C as a function of time: horizontal averages for ~1 period. In the different panels KF denotes results from 
testing the KuhfuB model and SW results from testing the Stellingwerf model. The results denoted by et,i include only the vertical 
velocity in the computation of the turbulent energy while for results denoted by et,2 the horizontal velocity was included as well. 


TDC 

convection zone 

overshooting region 


/9c 

C7 

/^max 

/^min 

Pc 

(7 

/^max 

/^min 

KF 

0.56 

23 

2.2 

0.07 

13 

1.6 

SW 

0.72 

25 

2.3 

0.22 

14 

1.7 


Table 2. /3 c -values for various convection models (/ $ Cj kf ~ a c,KF • OL r and ft Cj sw ~ a c,SW/ a r)- Peclet 
factors contain both the x- and ^-direction contributions to the turbulent velocity. For further details see 
text. 


tion from the vertical velocity component only. Our results 
indicate opposite behaviour in the convection zone and the 
overshooting region, respectively: in the convection zone, in¬ 
clusion of both components seems to lead to better results, 
whereas in the overshooting region inclusion of the vertical 
velocity component only seems to be preferable. That is true 
both for the KuhfuB and the Stellingwerf formulation. On 
the other hand, the inclusion of the Peclet correction fac¬ 
tor below the overshooting zone 77 seems to have a rather 
marginal influence on the quality of the results due to the 
small influence of the specific value of a r on the results, since 
a change in a r can effectively be compensated by a change 
in a c . This is an immediate consequence of, for example, 


Eq. [ 9 ] the relevant values of a r (see Table [TJ and the fact 
that in our case therefore a r • Pe C 1. 


3.1.4 A first analysis 

At first it might seem surprising that the inclusion of ra¬ 
diative losses does not significantly improve the predictions 
of the convective flux by the TDC models. The best fit nu¬ 
merical values for a c indeed change as a result of including 
a factor 77 ^ 1 , but there is hardly any improvement, if we 
compare the results for the convection zone with the over¬ 
shooting region. However, throughout most of the pulsation 
cycle the Peclet number Pe is of order unity or even less 
in the entire spatial domain analysed in our study. Indeed, 
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the horizontal average barely ever reaches the value of 0.75. 
Thus, we are dealing with convection in the low efficiency 
limit in both cases, i.e. for the convection zone proper and 
the overshooting region. For the He n zone most of the con¬ 
vective driving stems from the lowering of V a d rather than 
from the rather mild local minimum in K va a- Consequently, 
the correction due to radiative losses is similar in both re¬ 
gions and the correction factor 77 * cannot compensate the 
large difference found in the convective flux predictions for 
the convective and the overshooting region. Moreover, as it 
can be more easily seen from the form of ? 7 kf in Eq. 
it is the product a c a r only that enters into the calculation 
of the convective flux in the limit of small Peclet number. 
Thus, a change in a r just implies a scaling of a c throughout 
the simulation domain, since in the end we obtain a c from 
the ratio of Eq. § and Eq. ( |12| ). Hence, the problem of too 
large an overshooting found for the TDC models cannot be 
accounted for by just including a model of heat losses but 
must have an additional, more important contribution. - We 
also note that the inclusion of heat losses does at least pre¬ 
vent even more excessive overshooting underneath the He II 
partial ionisation zone which would otherwise be borne out 
by the Kuhfuss model. This follows from a comparison of 
Fig. |2(a)| and Fig. |2(b)| 

Possibly the disparity between the a c -faetors applicable 
to the convection or overshoot zone, respectively, we derive 
here would turn out to be even worse if similar work as 
the present one would be performed in 3D instead of 2D. 
In the 3D case overshooting is usually less pronounced (cf. 
Muthsam et al. 1995) than in 2D simulations with otherwise 


identical setup. We therefore expect an even larger ratio of 
the a c between the convection zone and the overshooting 
zone. This suggests that work in 3D should be undertaken 
in the future. 

We caution that the findings of the discussion above 
cannot be simply extrapolated to the convection zone caused 
by the partial ionisation of H 1 and He 1 . There, both the 
variation of K va a and the values of Pe are much larger. 

The inclusion of the horizontal component of the tur¬ 
bulent kinetic energy into the calculation of e t and, through 
Eq. § , into F c is more important, as is clearly demonstrated 
by the results in Table ^ The improvement within the con¬ 
vection zone for both TDC models is about 30%. But no such 
improvement is found for the overshooting zone: indeed, the 
model predictions become slightly worse. This can be un¬ 
derstood by means of the anisotropy of the flows: the ratio 
of horizontal to vertical kinetic energy is different for both 
regions, since vertical flow dominates inside the convection 
zone, while horizontal flow dominates outside. As the TDC 
models, implicitly or explicitly, assume some - usually con¬ 
stant - degree of anisotropy, this result demonstrates that 
this quantity is variable and in particular significantly differ¬ 
ent inside and outside the convection zone proper. Clearly, 
this is one route along which the current TDC models could 
be improved. 

It is interesting to approach the question of constancy of 
coefficients also from a different side for the He 11 convection 
zone of our model. It is a weak convection zone, and the 
Peclet numbers are small both in the convection zone proper 
and in the overshooting region. 

However, even if the anisotropy of the distribution of 
kinetic energy and the amount of overshooting could be 


taken care of, the problem of the variability of the qual¬ 
ity of the models with phase would likely remain, unless the 
convection models could incorporate the feedback between 
the mean structure and the turbulent kinetic energy as well 
as the convective flux in a more self-consistent manner. 

The fact that in the weak He 11 convection zone the 
parameters a c and a r enter only through 


Pc,. ■ = 


3>* • Pe* 


(13) 


as (3 c ,kf ~ ot c ,KF • OL r (KuhfuB) or /3 c ,sw ~ ol c ,sw / ot c ,sw 
(Stellingwerf), respectively, raises the question whether 
treating /3 C> * in the same way as has been done with <a c ,* 
(cf. Table 0 Fig.@ might not result in a better constancy 
of the flux predicted by the KuhfuB or the Stellingwerf pre¬ 
scription. Table [2] and Fig. [5] teach, however, that this is 
not the case. Rather, parameters such as ratio of the /3 Cj ,,- 
coefficients in the convection and the overshoot zone, respec¬ 
tively, or the variances tend to be worse than is true for the 
corresponding values based on the <a Cj *-coefficients. 

We turn attention to work of IGastine fe Dintransl 
(2011). Although they have performed their simulations for 
a case of idealised microphysics and plane parallel (as op¬ 
posed to spherical) geometry and differ, moreover, in the 
numerical setup, they have also found a temporal modula¬ 
tion of the convective flux (which can easily be read from 
our Fig. [TJ and, moreover, they have found that the convec¬ 
tive flux predicted by the Stellingwerf model is in somewhat 
better agreement with the simulations than those obtained 
from the KuhfuB model. The identified shortcomings of the 
most common TDC models used in the study of convection- 
pulsation coupling are thus quite likely of a generic nature 
and not specific to Cepheids of a very narrow part of the 
HR diagram. 

Concluding this section we want to consider the fol¬ 
lowing aspect. The o C 5 ma x/a: C5m in ratio refers to quantities 
integrated over the convection or overshoot zone. In no way 
could the predictions of the TDC models be interpreted in 
the pointwise sense (this has also not been suggested by the 
authors, as we hasten to remark). In this case the differences 
between the model predictions and the simulations would be 
even larger (by two orders of magnitude). From that point 
of view the averaging procedure still does a remarkable job 
in dealing with the point values which, taken per se, are 
practically meaningless. 


3.2 Energetics of the convection zone 

3.2.1 Basic quantities 


The energy equation of hydrodynamics (e.g. Castor (2004)) 
can easily be cast in the following form useful for our present 
purpose: 


Dt (e + ekin ) + (jp + £kin) div u — (14) 

— e div u — u • gradp + div (k grad T + ua) + pg • u 


or equivalently 

Dt + ~~) + (P + e kin) D t (t) = (15) 

— (ekin div u — u • gradp + div (k grad T + ua)) + g • u. 

P 
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Here, e denotes the internal and ekin the kinetic energy 
density, p is the (thermodynamic) pressure, p the density, 
a the stress tensor, k the conductivity, T the temperature 
and g the vector of gravitational acceleration. Dt is the La- 
grangian time derivative. Note that pDt(^) is close to the 
expression pdV of thermodynamics (instead of a change of 
volume there is now a rate of change, and we work with 
an energy density instead of an energy enclosed in a spe¬ 
cific volume, so that a term —edivtx appears on the right 
hand side to account for the rate of change of the volume) 
in Eq. ( fl4| . 

We now discuss several components of velocity and pres¬ 
sure, respectively, and their contribution to the so-called 
pdV-work. 

To begin with, we denote by /x the vector of momen¬ 
tum densities. In addition to the total velocity field u , we 
consider the fluctuating velocity field 

u , = /V_Ar ^\ ( 16 ) 

V p p p ) 

The bar denotes horizontal averaging, so that p r stands for 
the horizontally averaged r-component of the momentum 
density. Since the relevant quantities (/it, ix, p, p) are com¬ 
bined in a nonlinear way it is not possible to define the 
horizontal means in a fully consistent manner unless one is 
ready to deal with a multitude of emerging mixed terms. 
The detailed way of averaging has, however, practically no 
influence on the results in the present case due to the small¬ 
ness of the mixed terms. 

Various components of the total pressure are of physical 
interest, p denotes the thermodynamic pressure, as obtained 
from the equation of state. We then have 

P=P-P, (17) 

the fluctuating thermodynamic pressure. As to the macro¬ 
scopic pressure components, we consider 

Pt=^(Ur+ul) (18) 

and 

Pt = |((«r) 2 +«!)• (19) 


Pt is often referred to as the turbulent pressure whereas p t 
contains also the radial pulsational component. The term 
turbulent pressure is not meant to imply any statement 
about the degree of turbulence of the simulated flow in terms 
of the Reynolds number or similar quantities. We rather 
use it in the way in which it is used in the literature on 
one-dimensional simulations of stellar pulsation, denoting 
anything due to other than the radially averaged motions. 
Modelling the turbulent pressure is a key issue in ID recipes 


for the pulsation-convection interaction (e.g. Stellingwerf 
(1982), Yecko, Kollath, V Buchler (1998)). 

To be more precise consider that some of these models, 
for example those in the papers just cited, also include an 
eddy viscous pressure as more closely described in |Gehmeyr| 
(1992). This eddy viscous pressure is assumed to be due to 
turbulence generated by the shear of the radial pulsation; cf. 
the discussion around Eq. (4) and Eq. (5) in that paper. - 
In multidimensional simulations there is neither a need nor 
a really convincing way to distinguish between those two 
kinds of turbulent pressure. Consequently, p t refers to the 
total turbulent pressure in such models. 


3.2.2 The work components 

The Lagrangian form of the continuity equation, 


Dtp — — pdiv u. 


( 20 ) 


combined with Eq. (15) leads to work densities (per second) 
of the form pressure times divergence of velocity. The decom¬ 
positions of pressure and velocity (jp = p + p , u = u + u') 
results in a number of work densities. Some of these are 
significant for the present discussion as inspection of data 
shows. 

The cycle averages over the horizontal average p • div u 
of these work densities ( Buchler V Kolldth| pOOO)) are the 
work integrands and express the amount of energy con¬ 
verted into pulsational energy. In Fig.[6]the work integrands 
are plotted for the fluctuating velocity. When comparing to 
Fig-H it can be seen that for the thermodynamic pressure 
dampening occurs at the top of the convection zone and that 
the transition from excitation to dampening is within the 
overshooting zone. The turbulent pressure has a dampening 
effect throughout where it is not small. 

These work densities can then be integrated, at any 
time, over the convection zone (c), the overshooting zone 
(o). For the present purpose the definition of the convection 
zone is the same as given in Sect. |3.1.2| The overshooting 
zone is now defined a little bit differently from the one pro¬ 
vided there. Namely, now there is no gap between it and the 
convection zone: points below the convection zone belong to 
the overshooting region if they satisfy T 3 < 1.51 and fulfill, 


V — V a d < 0 (instead of V — V a d < —0.03 in Sect. 3.1.2 ) 
This results in terms such as, for example, 


W°(p,u') = 


') = / p{x) div u'(x)dx, 
J o 


referring to the overshooting zone. In this way the rather 
different physical nature of the zones is taken into account. 

Within such an integral positive and negative contribu¬ 
tions may cancel out very significantly. This can be seen by 
considering the integrated modulus of the pdV-terms such 
as 



| p(x) div u {x)\dx 


alongside with the original values. 

For the thermodynamic pressure p the variation with 
phase is shown in Fig. [7] The work W (p,u) can be seen 
to mainly follow the pulsation pattern for each zone, more 
clearly so in the overshooting zone. The average values 
W (p, it) (see Tab. ^ are positive everywhere, since exci¬ 
tation work can be observed at the top of the overshooting 
zone (see Fig. [7] (a)). In the convection zone the pdV-terms 
reach their smallest values during contraction while the av¬ 
erage values are always positive and the modulus shows little 
variation indicating roughly constant activity. The panels on 
the right-hand side of Fig. [7] show that in the overshooting 
zone the modulus of the pdV-terms derived from fluctuation 
quantities reach their highest values right after maximal con¬ 
traction at the beginning of the expansion. Cancellation ef¬ 
fects in the work integral are most dramatic for the W(p, u') 
terms, reaching up to about three orders of magnitude (com¬ 
pare Fig[7](c) and (d)). 

For the turbulent pressure the pdV-terms follow the 
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250 300 350 400 450 



Figure 6. The work integrands for the fluctuating velocity field, 
included, while in the right panel only the pressure perturbations 


250 300 350 400 450 



(b) 


the left panel the pulsational component of each pressure term is 
used. The top of the simulation box is to the left of each panel. 


region thermodynamic pressure turbulent pressure 



W (p, u) 

W(p, u’) 

W(p',u) 

W(p',u') 

W (p t , u) 

W(p t ,u') 

W(p' t ,u) 

W(p' t ,u') 

convection zone 
overshooting region 

3.7e+27 

9.3e+27 

2.6e+27 

-5.0e+26 

7.4e+25 

5.8e+26 

4.0e+25 

4.5e+26 

-2.9e+27 

-2.5e+27 

-7.5e+26 

-2.3e+27 

-8.2e+26 

-2.9e+27 

-5.1e+26 

-2.4e+27 


Table 3. Average values of the components of the work integral in various regions. Note that for the turbulent pressure all 
values are negative. 


pulsation pattern if at least one component includes pulsa¬ 
tion. In the case of W (p' t ,u f ) one sees (Fig. U (§) and ( h )) 
that the dampening effects of the turbulent pressure per¬ 
turbation are at full contraction at their maximum in the 
overshooting zone but at their lowest in the convection zone. 
The modulus W ^ (p' t ,u') in the convection varies little. 


4 DISCUSSION AND CONCLUSIONS 


While simulations of pulsating stars (Cepheids, RR Lyr) 
have been computationally one-dimensional for decades, as¬ 
suming radial symmetry, multidimensional calculations start 
to be feasible only recently. They allow in particular the 
investigation of the pulsation-convection interaction which 
is considered a main open problem in the simulation of 


such stars (Buchler (2009)). We have constructed a two- 
dimensional model regarding convection and pulsation in a 
Cepheid type star. The justification for adopting two di¬ 
mensions at the present state of the matter can be taken 
from a comparison of 2D and 3D results in Gero ux &; De-| 
upree| ( 2014| . These results show that (for models off RR 
Lyr stars) many quantities such as peak convective flux as 
a function of phase, growth rates or amplitudes are quite 
similar in 2D and 3D, differences being mainly visible in the 
slope of the light curve at various times in the descending 
branch. The aim of our investigation is to figure out basic 
properties of convection in the presence of pulsation and to 
draw inferences on existing ways of modelling these issues 
via time dependent convection models (TDCs) in a compu¬ 
tationally ID setting. 


4.1 Convection modelling in pulsating stars 

The TDCs just mentioned entail a number of coefficients 
not priorly known. With respect to the coefficient a c , the 
calibration factor for obtaining the physically correct con¬ 
vective flux values from the TDC models of Kuhfufi and 
Stellingwerf, the following conclusions can be drawn from 
our work. 

In the He II partial ionisation zone and its surrounding 
overshooting zone convection is subject to large radiative 
losses, as is witnessed by a low Peclet number (Pe < 1 pre¬ 
dominantly). It is thus difficult to calibrate a c separately 
from a coefficient a r , connected to radiative losses, at least 
from our present model. After all, as explained in Sect. 3.1, 
these coefficients enter essentially only in the combination 
a c ■ a r (or a c /a r in case of the Stellingwerf model) for our 
case of low Peclet numbers. 

Our models show that using one and the same coef¬ 
ficient a c for the convection zone (c.z.) and the overshoot 
zone (o.z.) would result in much too large overshooting pre¬ 
dicted by the TDCs. However, a c and a r enter only in the 
combination described in Sec. |3.1.2] Consequently, including 
radiative losses by any choice of a r in the prescription for the 
convective flux cannot correct the discrepancy of too large 
fluxes predicted by the model expressions for the overshoot¬ 
ing zone in comparison with the interior of the convection 
zone. At present, the main effect of accounting for radiative 
losses, i.e. introducing the factor p in Eq. (3), is to at least 
limit a more excessive overshooting which would otherwise 
result from applying the prescription of Kuhfufl (see Fig. 2). 
This model is more affected by that problem than the model 
of Stellingwerf. 

Accounting for the change of anisotropy of the contri- 
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Figure 7. The work integrals including the thermodynamic pressure for ~ 1 period. The components used are indicated on the vertical 
axis of the different panels. 
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Figure 8 . The work integrals including the turbulent pressure for ~ 1 period. The components used are indicated on the vertical axis 
of the different panels. 
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butions to turbulent kinetic energy, however, is much more 
important. To avoid excessive overshooting in the model re¬ 
quires introducing a functional dependence in Eq. (3) which 
can take care of the fact that convective zones show stronger 
vertical fluxes than horizontal ones while the opposite holds 
for overshooting zones. 

However, as we have shown, even after such corrections 
a large discrepancy remains. This requires to at least vary 
a c5 * by a factor of 1.5 for the Stellingwerf model and a factor 
of 4 for the KuhfuB model to obtain convective fluxes from 
the models which are consistent with our numerical simula¬ 
tions. In addition, there is a clear systematic trend of a c as 
a function of depth within the convection zone, see Fig. [3] 
The effect of this trend on modelling the pulsation in the 
traditional way is unknown. 

Some contribution to this might actually be due to the 
coupling between the mean structure and the turbulent ki¬ 
netic energy which is certainly more complex than expressed 
through Eq. (3). 

We should mention that the present comparison is al¬ 
ready the most favorable one for the TDC models in the 
sense that we only take the model expressions to be meant 
to predict statistical averages: variations for particular in¬ 
stants of phase for a specific point in time may be larger 
and any attempt to interpret the expressions in a point-wise 
sense is falsified immediately by unacceptably large discrep¬ 
ancies (at least two orders of magnitudes larger than the 
already large differences observed for the averaged quanti¬ 
ties considered here). 

In this sense the present simulations unveil major weak¬ 
nesses of TDC models used to explain the pulsational prop¬ 
erties of Cepheids. The present simulations allow to quantify 
the discrepancies expected from the TDC models for a cer¬ 
tain group of Cepheids. Since idealized simulations such as 


those of Gastine &; Dintrans (2011) come to similar con¬ 


clusions on a qualitative level, the discrepancies found are 
most likely of a fundamental nature and thus have to be 
taken care of. 

Otherwise one may fall into the trap of obtaining the 
right results for what might be the wrong reasons. Our re¬ 
sults suggest that in ID simulations one gets artificially large 
overshooting below the convection zone when applying iden¬ 
tical values of a c in the c.z. and the o.z.. This reinforces the 
doubts which have been expressed by |Smolec &; Moskalik| 
(2008) in respect to the achievements on double mode pul¬ 


sation in ID Cepheid models. These doubts are based on 
the suspicion of too large overshooting in the models which 
have shown double mode pulsation in a broader parameter 
range than previously. 

One additional point is that 2D simulations tend to lead 
to larger overshooting than similar 3D simulations produce 
(Muthsam et al. (1995)). Considering this aspect as well, the 
problem of depth-dependency of a c may well be even more 
severe than is borne out by the present investigation. 


4.2 Energetics of the convection zone 

Our results regarding the contribution of fluctuating quanti¬ 
ties to the work integral reveal an intricate change of corre¬ 
lations as a function of phase. In Table [3] we compare values 
of the contributions to the work integral, integrated over the 
convection zone, and the overshoot zone. Considering contri¬ 


butions to the work integral where both entering quantities 
are fluctuating ones (instead of corresponding to the radial 
symmetric component) it turns out that no such combina¬ 
tion can be disregarded: at least for one domain of integra¬ 
tion every combination leads to a contribution of the same 
size as the others. That means that in principle each TDC 
model should faithfully represent the terms corresponding 
to all of the possible combinations (plus those with just one 
fluctuating quantity of course). It can be expected that work 
integral contributions with two fluctuating quantities will be 
more important in cases of stronger convection zones, since 
due to the more pronounced nonlinearity, their higher order 
terms will likely be larger. 

A comparison of our simulations to TDC models re¬ 
garding the energetics of the pulsation-convection interac¬ 
tion cannot be undertaken in the same manner as has been 
possible for the o-factors discussed above. There, the rele¬ 
vant flux as delivered by the TDC models can, by a sim¬ 
ple algebraic formula, directly be evaluated and the scaling 
factor be determined; the constancy or variability of this 
factor then immediately leads to an assessment of the qual¬ 
ity. This is to be contrasted to the case of the work integral. 
There, the relevant quantities in the TDC (turbulent energy 
or pressure plus furthermore others, not really provided by 
the usual TDCs) are not algebraic in nature. Instead, the 
equation for turbulent energy is an evolution equation where 
et attains a value deemed physically correct only after it 
has been evolving for some time. Hence, comparison to the 
quantities from the multidimensional simulations cannot be 
performed on the basis of a temporal snapshot. At the same 
time it is not really possible to include an equation for the 
turbulent energy et in the multidimensional simulations: the 
simulations explicitly calculate et while the equation for tur¬ 
bulent energy in the TDCs models et . Consequently, there is 
a branching quite at the roots of the modelling processes be¬ 
tween one- and multidimensional simulations which prevents 
a meaningful comparison of the two values of et by solving 
the TDC equation alongside in multidimensional work. Af¬ 
ter all, there is no room to allow for backreaction of the 
TDC variant of et in multidimensions. - Our calculations 
however show that, for modelling these values in the sense of 
TDCs, phase dependent properties of hydrodynamic quanti¬ 
ties ought to be built into the models. After all, the contribu¬ 
tions to the work integral (which are related to correlations 
of physical quantities) have a distinctive phase dependence 
which, moreover, is not in all cases similar in the convection 
and the overshoot zone. This cannot really come as a sur¬ 
prise given that convection in a pulsating star is not in a 
statistically steady state. Rather, it continually evolves in a 
specific way: at first convection cells form at the top of the 
convection zone at a certain phase, plumes then descend, 
overshoot into the stable zone below at some later phase 
and ultimately the remnants die off deep below the convec¬ 
tion zone. Explicitly or implicitly, each TDC should account 
for such a sequence of events. 


4.3 Requirements for time dependent convection 
models 

The simulations reported above suggest some guidelines for 
the development of revised TDCs to eventually be used in 
ID simulations. Radiative losses have to be included in the 
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model, but in addition the convective flux will have to evolve 
in an explicitly phase dependent manner. Therefore an ad¬ 
ditional evolution equation for F c seems inevitable, similar 


to the model of Gough (1977), further developed in Balm- 


forth (1992), or the rarely used 3-equation variant of the 


KuhfuB model. Moreover, the distribution of kinetic energy 
clearly changes as a function of depth in a non-trivial man¬ 
ner and most likely also features phase dependencies. Thus, 
a fully non-local Reynolds stress model as given by |Canuto| 
(1993), Canuto (1997) and Canuto (1999) may be required 
as a minimum level of complexity to substantially improve 
TDC models in comparison to those presently used. 


4.4 Further avenues of research 


When finally turning our attention to further work along 
these lines, several items come to mind. One obvious goal 
consists in the construction of Cepheid models with a 
stronger He n convection zone in order to more fully sam¬ 
ple the parameter space. Such material may be crucial for 
improved TDCs for use by the ID pulsation modelling com¬ 
munity. 

One further goal is an investigation of the upper, 
H+Hel, convection zone and, hence, the atmosphere. We 
have already performed some initial runs with the neces¬ 
sary resolution which is considerably higher in the upper 
layers than in the models which we have reported here. Due 
to the fact that, with higher resolution, the H+Hel convec¬ 
tion zone is much stronger than in our present models we 
have had some difficulty with the higher resolution models. 
In particular, throughout the atmosphere (and also below) 
there is vivid non-radial shock activity. This has lead to nu¬ 
merical problems which required the development of numer¬ 
ically suitable boundary conditions. Such conditions were 


devised in the context of solar granulation ( Grimm-Strele et 
al. (2015)). A modification of these conditions has recently 
been found to work for the Cepheid case. Shock activity 
with shocks up to a Mach number of about 4, triggered by 
convection, has also been found in preliminary calculations 
we have done with these new conditions. This highly dy¬ 
namical situation in the atmosphere of the Cepheid model 
also leads to a highly corrugated optical surface. To describe 
such a situation by a static or even a time-dependent model 
one-dimensional seems extremely challenging, and the impli¬ 
cations for an analysis of photometric or spectroscopic data 
will need to be figured out. 

A different avenue of research could center on investing 
pulsation in a ring (or, in 3D, a spherical shell) as initiated 
in preliminary work ( |Kup ka et al. ( 2014| ). Such simulations 
could shed light on the excitation of non-radial modes if 
they indeed show up, connected or not connected with con¬ 
vection, whatever the outcome of the simulations would be. 
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